/*!
 **************************************************************************
 * \file convdiff1.cpp
 * This code simulate Natural Convection in a unit square.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz [ Sat May 23 12:06:36 CDT 2009 ]
 ***************************************************************************
 * The equations to be solved are of the form:
 * \f[
 * \frac{\partial \phi}{\partial t} + \frac{\partial}{\partial x_j} 
 * \Big(u_j\phi\Big) =
 * \frac{\partial}{\partial x_j}\left(\Gamma\frac{\partial\phi}{\partial x_j}\right) + S,
 * \,\,\,\,\,\,\,\,\textrm{for  } j = 1,2,
 * \f]
 *
 * where \f$ \phi \f$ is a scalar variable \f$T, \rho, u_1\f$ or \f$u_2 \f$.
 * \f$ \Gamma \f$ is a diffusion coefficient and \f$ S \f$ is the source term.
 *
 * In this case we have an energy equation in terms of temperature \f$T\f$, 
 * coupled with the Navier-Stokes equations. 
 * The coupled equations are solved with the SIMPLEC strategy.
 * 
 * These equation are solved in a unit cube 
 * \f$ x, y \in [0,1] \times [0,1] \f$ as shown in
 * the next figure:
 * \image html convection.png "Domain in 2D" width=10cm
 * \image latex convection.eps "Domain in 2D" width=10cm
 * Two opposite walls have constant temperatures but different. The other walls
 * are adiabatic. The no-slip condition is applied on all walls.
 *
 *  To compile and run this example type the next commands:
 * \par
   \verbatim
   % make
   % ./convdiff1 \endverbatim
 * 
 ***************************************************************************
 */


// Note: Schemes: Quick_..., CDS_..., Upwind_...
#include <iostream>
#include <cmath>

#include "Meshes/Uniform.hpp"
#include "Equations/Momentum_XCoDi.hpp"
#include "./LDS/Upwind_XLD.hpp"
#include "Equations/Momentum_YCoDi.hpp"
#include "./LDS/Upwind_YLD.hpp"
#include "Equations/ScalarEquation.hpp"
#include "./LDS/CDS_LD.hpp"
#include "Equations/PressureCorrection.hpp"
#include "Schemes/Simplec.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/diff_operators.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;

// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

template<class T_x, class T_y, class T_p>
double SIMPLEC(T_x &, T_y &, T_p &, int, double, double, double);

double errorT, errorX, errorY;

void inlet_outlet(ScalarField2D &u, ScalarField2D &s, double dx, double dy) 
{
  double y;
  int top = u.ubound(secondDim);
  int mid_y = top / 2;

  // INLET
  for(int j = mid_y; j < top; ++j) {
    y = j * dy - 0.5;
    u(0, j) = 24 * y * (0.5 - y);
    s(0, j) = 2 * y * y * (3 - 4 * y);
  }  

  // OUTLET
  for(int j = 0; j < top; ++j) {
    y = j * dy - 0.5;
    u(u.ubound(firstDim), j) = 0.75 * (1 - 4 * y * y);
    s(s.ubound(firstDim), j) = 0.75 * y  - y * y * y + 0.25; 
  }
}

int main() 
{

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    int stop;
    plotter("set grid");
    plotter("set contour");
    plotter("set cntrparam level discrete 0.4, 0.3, 0.2, 0.1, 0, -0.1, -0.2, -0.3, -0.4");
    plotter("set view 0,0");
    plotter("set size 0.73,1");
    plotter("unset surface");
#endif

    double length_x, length_y;
    int num_nodes_x, num_nodes_y;
    double tolerance, dt;    
    int max_iter, max_time_steps, frequency;
    double Reynolds, Gamma;

// Reading data from the input file
    std::ifstream input_cin ("input1");
    input_cin >> length_x
	      >> length_y
	      >> num_nodes_x 
	      >> num_nodes_y 
	      >> dt
	      >> max_time_steps
	      >> tolerance
	      >> max_iter
	      >> Reynolds
	      >> frequency;
    input_cin.close();

    Gamma = 1.0 / Reynolds;

//
//  Mesh definition
//
    StructuredMesh<Uniform< double, 2> > mesh(length_x, num_nodes_x,
					      length_y, num_nodes_y);
    mesh.print();
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);
    int num_vols_x = mesh.getExtentVolumes(X);
    int num_vols_y = mesh.getExtentVolumes(Y);

    ScalarField2D p(mesh.getExtentVolumes());      // pressure
    ScalarField2D pp(mesh.getExtentVolumes());     // pressure correction
    ScalarField2D u(mesh.getExtentVolumes());      // u-velocity
    ScalarField2D v(mesh.getExtentVolumes());      // v-velocity
    ScalarField2D us(num_nodes_x, num_vols_y); // u-velocity staggered
    ScalarField2D vs(num_vols_x, num_nodes_y); // v-velocity staggered

    ScalarField2D pn(mesh.getExtentNodes());
    ScalarField2D un(mesh.getExtentNodes());
    ScalarField2D vn(mesh.getExtentNodes());    
    ScalarField2D vort(mesh.getExtentNodes()); 
    ScalarField2D streamn (mesh.getExtentNodes());
    
// Boundary conditions
    inlet_outlet(un, streamn, dx, dy);
    NumUtils::staggerFromNodesX(us, un);

#ifdef WITH_DX
    InOut::writeToFile_DX(pn, 0, "./Data1/pres.", dx, dy);
    InOut::writeToFile_DX(streamn, 0, "./Data1/stre.", dx, dy);
    InOut::writeToFile_DX(un, vn, 0, "./Data1/velc.", dx, dy);
#endif

// Linear System
    SparseMatrix< Diagonal<double, 2> > A(num_nodes_x, num_nodes_y); 
    ScalarField2D             b(num_nodes_x, num_nodes_y);

//
// Xmomentum equation definition
//
    Momentum_XCoDi<Upwind_XLD<double, 2> > mom_x(us, A, b, mesh.getDeltas());
    mom_x.setGamma(Gamma);
    mom_x.setDeltaTime(dt);
    mom_x.setDirichlet(LEFT_WALL);
    mom_x.setDirichlet(RIGHT_WALL);
    mom_x.setDirichlet(TOP_WALL);
    mom_x.setDirichlet(BOTTOM_WALL);
    mom_x.setUvelocity(us);
    mom_x.setVvelocity(vs);
    mom_x.setPressure(p);
    mom_x.print();

//
// Ymomentum equation definition
//    
    Momentum_YCoDi<Upwind_YLD<double, 2> > mom_y(vs, A, b, mesh.getDeltas());
    mom_y.setGamma(Gamma);
    mom_y.setDeltaTime(dt);
    mom_y.setDirichlet(LEFT_WALL);
    mom_y.setDirichlet(RIGHT_WALL);
    mom_y.setDirichlet(TOP_WALL);
    mom_y.setDirichlet(BOTTOM_WALL);
    mom_y.setUvelocity(us);
    mom_y.setVvelocity(vs);
    mom_y.setPressure(p);
    mom_y.print();
//
// Pressure correction equation definition
// 
    PressureCorrection<Simplec<double, 2> > press(pp, A, b, mesh.getDeltas());
    press.applyBounds(1, num_nodes_x-1, 1, num_nodes_y-1);
    press.setUvelocity(us);
    press.setVvelocity(vs);
    press.set_du(mom_x.get_du());
    press.set_dv(mom_y.get_dv());
    press.setPressure(p);
    press.print();

    ofstream sorsumfile("./Data1/sorsum"),
      errorXfile("./Data1/errorX"), errorYfile("./Data1/errorY");
    
    double sorsum;
    int iteration = 0;

    ScalarField2D stream (mesh.getExtentVolumes());
    ScalarField2D vortv (mesh.getExtentVolumes());
    NumUtils::interpolateToVolumes(stream, streamn);

    ScalarEquation< CDS_LD<double,2> > laplace(stream, A, b, mesh.getDeltas());
    laplace.setGamma(1);
    laplace.setNeumann(TOP_WALL, 1);
    laplace.setDirichlet(BOTTOM_WALL);
    laplace.setDirichlet(LEFT_WALL);
    laplace.setDirichlet(RIGHT_WALL);
    laplace.setUvelocity(vortv);

    //    NumUtils::interpolateToNodes(pn, p);
    //    NumUtils::interpolateToNodesUs(un, us);
    //    NumUtils::interpolateToNodesVs(vn, vs);
  
    for(iteration = 1; iteration <= max_time_steps; ++iteration) {
      sorsum = SIMPLEC(mom_x, mom_y, press,
		       max_iter, tolerance, dx, dy);

      if ( (iteration % frequency) == 0 ) {
	NumUtils::interpolateToNodes(pn, p);
	NumUtils::interpolateToNodesUs(un, us);
	NumUtils::interpolateToNodesVs(vn, vs);	    
	DiffOperators::vorticity(vort, un, vn, dx, dy);

	NumUtils::interpolateToVolumes(vortv, vort);

	laplace.calcCoefficients();	
	Solver::TDMA2DY(laplace, tolerance, max_iter);
	laplace.update();

	NumUtils::interpolateToNodes(streamn, stream);	    

#ifdef WITH_DX
	InOut::writeToFile_DX(pn, iteration, "./Data1/pres.", dx, dy); 
	InOut::writeToFile_DX(un, vn, iteration, "./Data1/velc.", dx, dy);
	InOut::writeToFile_DX(vort, iteration, "./Data1/vort.", dx, dy);
	InOut::writeToFile_DX(streamn, iteration, "./Data1/stre.", dx, dy); 
#endif

#ifdef WITH_GNUPLOT
	InOut::writeToFile(pn, 1, "./Data1/pres.", dx, dy); 
	InOut::writeToFile(un, 1, "./Data1/uvel.", dx, dy);
	InOut::writeToFile(vn, 1, "./Data1/vvel.", dx, dy);
	std::cout << "\n >---> Press <1> and then <enter> to continue \n";
	std::cin >> stop;	
	//	    plotter("splot \"Data1/pres.00001\" w l");
	plotter("splot \"Data1/uvel.00001\" w l");
	//	    plotter("splot \"Data1/vvel.00001\" w l");
#endif
      }
	
      sorsumfile << iteration << "\t" << sorsum << endl;     
      errorXfile << iteration << "\t" << errorX << endl;
      errorYfile << iteration << "\t" << errorY << endl;
      cout << "=> Iter = " << iteration 
	   << " e_X = " << errorX
	   << " e_Y = " << errorY
	   << endl;
    }

#ifdef WITH_DX
    std::cout << "\n\n NOTE: The data were written in the [./Data1] directory in DX format";
#endif

    return 0;
}


template<class T_x, class T_y, class T_p>
double SIMPLEC(T_x &mom_x, T_y &mom_y, T_p &pressure, 
	       int max_iter, double tolerance, double dx, double dy) 
{
    double sorsum = 10.0, tol_simplec = 1e-2;
    int counter = 0;

    while ( (sorsum > tol_simplec) && (counter < 20) ) {
	mom_x.calcCoefficients();	
	Solver::TDMA2DY(mom_x, tolerance, max_iter);
	errorX = mom_x.calcErrorL1();
	mom_x.update();

	mom_y.calcCoefficients();
	Solver::TDMA2DY(mom_y, tolerance, max_iter);
	errorY = mom_y.calcErrorL1();
	mom_y.update();

	pressure.calcCoefficients();
	Solver::TDMA2DY(pressure, tolerance, max_iter);
	pressure.correction();
	sorsum = fabs( pressure.calcSorsum() );

	++counter;
    }

    return sorsum;
    
}
